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Within a fully microscopic setting, we derive a variational principle for the non-equilibrium steady 
states of chemical reaction networks, valid for time-scales over which chemical potentials can be 
taken to be slowly varying: at stationarity the system minimizes a global function of the reaction 
fluxes with the form of a Hopfield Hamiltonian with Hebbian couplings, that is explicitly seen to 
correspond to the rate of decay of entropy production over time. Guided by this analogy, we show 
that reaction networks can be formally re-cast as systems of interacting reactions that optimize the 
use of the available compounds by competing for substrates, akin to agents competing for a limited 
resource in an optimal allocation problem. As an illustration, we analyze the scenario that emerges 
in two simple cases: that of toy (random) reaction networks and that of a metabolic network model 
of the human red blood cell. 



INTRODUCTION 

The dynamics and thermodynamics of chemical reac- 
tion networks is a subject that goes back at least to [THE]. 
Recent years have seen considerable interest in the prob- 
lem at different levels, from the characterization of their 
mass-action kinetics (7] and of their stochastic thermo- 
dynamics based on the chemical master equation (8] [9] , 
to the analysis of their non-equilibrium steady states 
(NESS) [TUlll3| . Besides their relevance for fundamen- 
tal understanding, these approaches provide an impor- 
tant frame for the study of real biochemical systems, like 
genome-scale reconstructions of cellular metabolic net- 
works [HHTH]. The most basic information about these 
systems is usually encoded in the matrix of stoichiometric 
coefficients, representing in essence the (weighted) topol- 
ogy of the couplings between chemical species and reac- 
tions. With uncertainties about kinetic parameters and 
reaction or transport mechanisms often preventing large- 
scale kinetic approaches (with exceptions like the the 
metabolism of human erythrocytes [T7]), the challenge 
at the simplest level is that of building stoichiometry- 
based predictive models of metabolic activity at genome 
scale. Much information on the organization of reaction 
fluxes in NESS can indeed be obtained from constraint- 
based models that rely on minimal mass-balance [T8H20] 
or stability [2TJ [22] assumptions. Such descriptions re- 
volve around pre-defined (and physically motivated) sets 
of local constraints, enforcing for instance mass balance 
at each metabolite node in the network. It would be 
interesting (and instructive) to derive the relevant local 
constraints as the result of a mathematical analysis of 
the dynamics taking place on the network, both to fur- 
ther clarify the assumptions behind the models and to 
highlight their limitations. 

This work revisits the joint dynamics of concentrations 
and reaction fluxes in chemical networks in a statistical 
mechanics frame. In specific, we obtain a variational 



principle that relates fluxes (i.e. the average number 
of microscopic transitions per time per volume for each 
process) in NESS to the minima of a global function 
where stoichiometric coefficients and steady-state con- 
centrations appear as parameters. This function is remi- 
niscent of the Hamiltonian of a Hopfield model with Heb- 
bian couplings [23] , with stoichiometry playing the role of 
the patterns. An analysis of its physical meaning explic- 
itly shows that reaction networks dynamically converge 
towards states where the use of the available compounds 
is optimized and the rate of decay of entropy production 
is minimized. The flux organization problem turns out to 
have remarkable similarities with that of optimal resource 
allocation by heterogeneous agents, as described e.g. by 
Minority Games [24l [25] . Systems of this type generically 
undergo a transition from an ergodic phase (the NESS is 
independent of the initial conditions of the dynamics) to 
a non-ergodic one when the ratio between the number of 
reactions and that of chemical species is changed. In our 
case, the two regimes are described by different sets of 
local constraints. We shall first explore this scenario in 
toy "random" reaction networks where such a transition 
can be fully analyzed numerically. Then a simple real 
system will be considered, namely the reduced metabolic 
network model of human erythrocytes. Finally, we shall 
discuss the relevance of these results for the quantitative 
analysis of cellular metabolism. 

ANALYSIS 

The variational principle 

We consider an open system enclosed in a volume V 
(a reactor or "cell", for brevity), formed by M distinct 
chemical species (labelled fx) that can be processed by 
N distinct reactions (labelled at fixed temper- 

ature T, pressure, ionic strength and pH. The reaction 
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stoichiometry is described by the coefficients with 
the convention that negative (resp. positive) coefficients 
identify substrates (resp. products) in the 'forward' di- 
rection of reaction i. Individual reaction events occur 
stochastically with rates (number of events per unit time) 
proportional to the substrate concentrations, which vary 
in time. At stationarity and for ideal systems, the Gibbs 
energy (GE) change per mole associated to each reaction 



fe M > 0) between time t and time t + St. Taking V to be 
fixed, the time-evolution of concentrations is simply 



x»(t + 5t) -x^(t) = 5x^{t) = 
We now focus on the quantity 



5n^(t) 
V 



y i (t) = y i , -J2^log{x fi {t)/x } 



(4) 



(5) 



AG, = AG 2 , + RTj2ti log(*7zo) 



(1) 



(where R is the gas constant, x M the intracellular con- 
centration of species ft,xga reference concentration and 
AG^o the GE change in standard conditions at con- 
centration xq) characterizes its distance from detailed 
balance. Specifically, the forward {4>i,+) and reverse 
(<fii -) fluxes of i, i.e. the average number of transitions 
per time per volume, satisfies the relation —/3AGi = 
log(0i 1+ /0i _) with 1//3 = RT [10]. At equilibrium, 
<f>i t+ = and AG; = for each i. 

We focus on the time evolution of the internal composi- 
tion of such a chemical reactor. The dynamics of this sys- 
tem is driven essentially by two factors: (a) the fact that 
molecules can stochastically enter or leave the system, 
and (b) the fact that reactions events occur (stochasti- 
cally) inside the system. Reasonably, then, the NESS will 
depend (a) on the rate at which molecules can cross the 
system's boundaries per unit volume (intake or outtake 
fluxes, denoted by m m ), and (b) on the rates of the internal 
reactions (more precisely, on the average net number of 
microscopic transitions per time per volume, denoted by 
4>i). In the following we will show specifically that, given 
the boundary fluxes (characterizing the environment), 
over time scales for which the chemical potentials can be 
assumed to vary slowly (so that concentrations can be 
assumed to be roughly constant) the internal fluxes in a 
NESS minimize the function 



> 



(2) 



where x^ denotes the concentration of species fi. 

Assume that at time t = the system is characterized 
by molecular populations (number of molecules) n M (0). 
Consider a time interval of size St and let Vi (t) denote the 
net number of transitions of reaction i that take place be- 
tween time t and time t+St. The latter is a random num- 
ber governed by a probability law (which shall depend, 
e.g., on the concentration of substrates) that we leave un- 
specified for the moment. The corresponding variation of 
7i A "s is given by 



n M (t + St) - n"(t) = Sn^it) 



lesw -&"(*) , (3) 



where & M (i) stands for the net (random) number of 
molecules of species fi taken in (if < 0) or out (if 



(with j/i.o = — jSAG^o), whose change between times t 
and t + St is given by 

Vi(t+5t)- yi (t) = 5vi(t) = - log . (6) 

n 

with initial conditions yi(0) = yi,o — J2fj. £f log^W/zo]- 
Equation ^ contains sources of stochasticity in the Wf's 
and the 6 M 's. As a consequence, the "macroscopic" vari- 
ables x^ and yi will fluctuate stochastically as well. Our 
aim is to characterize the steady state(s) of We make 
the following simplifying assumptions: 

Al. Molecular populations are large enough to allow 
us to treat x^ as a continuous variable. This is 
generically assumed to be the case in real cells, 
although e.g. in E. coli the number of copies of 
certain small molecules can be as low as a few tens 
(corresponding to a concentration of the order of 10 
nM [37]). The effects induced by molecular noise 
can be non trivial [9] and accounting for it might 
alter the emerging picture [55] . 

A2. The quantity Sx^/x^ ~ x^St/x^ is small (i.e. the 
chemical potential g^ = g$ +RT log x M , with gft the 
standard chemical potential, is roughly constant), 
so that Syi(t) can also be taken to be small. 

Under the above assumptions the right-hand side of ^ 
is easily linearized to yield 



Vi{t+8t)-yi(t) 



V ^ 

j 



v 



Si Sj 



-E 

V ^ 



x**(t) 
(7) 



This equation highlights the way in which concentrations 
affect the time evolution of the system and, in princi- 
ple, one would now need to analyze the coupled system 
formed by Q and However, for simplicity, we re- 
place x^it) with some time-independent limit value x^. 
This approximation can only be justified as long as one 
considers evolution over time scales shorter than x^/x 11 . 
If however concentration changes are sufficiently slow (in 
agreement with homeostasis) it is reasonable to expect 
that it will hold over time scales much longer than those 
required to reach a NESS. With this, Q takes the form 



Vi{t + 6t)-yi(t) 
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where the "couplings" Jij are defined as 



Jij — 



T P 



(9) 



The steady state of ^ can now be obtained straightfor- 
wardly by dividing both sides by St and averaging over 
time. This gives 



5t 



7 / 



V5t 



er / b»gt) 

VSt 



(10) 



Note that (vi(t)/VSt) = fa is the net flux of reaction i 
(the average net number of microscopic transitions per 
time per volume) while (t) (t) / V^<5i) = u M represents 
the uptake of species ii (the average number of molecules 
of species \x per time per volume entering or leaving the 
cell). Defining 



we therefore have 

Syi 
6t 



1 dH 



(11) 



(12) 



with 



-1 2 



> 



(13) 



This implies that, on average, the stochastic dynamics of 
reactions collectively minimizes the function H so that 
the NESS of ( 12 ) correspond to the solutions of the min- 



imization problem 



min H 



(14) 



constrained by the uptake values and by the fact that the 
fa's are assumed to be bounded by enzyme availability, 
so that fa 



J i. ni/n , .max \ 



Following [241 G2] , we can characterize the behavior of 
fluxes in NESS by noting that the solutions of (12 1 for 
t — » oo (which are expected to be linked to the Gibbs 
energy change of reaction i in a NESS) are generically of 
the form (yi) ~ 7^ with 7, = hi + = Jij<ftji since ji is 
constant in NESS. In practice, the stochastic dynamics of 
Ui for different reactions will be characterized by different 
values of 7^, depending on the asymptotic behavior. If yi 
tends to a finite value as t — > 00 then 7^ = 0. Recalling 
that the Gibbs energy change is related to the forward-to- 
reverse flux ratio through the detailed balance condition, 
this describes the case of a reaction whose microscopic 
transitions occur bidirectionally even for t — > 00, i.e. such 
that the forward and reverse fluxes are both non-zero in 
the steady state. Specifically, its flux is determined by 



the condition hi = — J2j Jij<t>j- When 7i 7^ 0, instead, j/j 
increases or decreases linearly in time (after a transient) 
and diverges for t — > 00 so that, in the corresponding 
NESS, reaction i will be characterized by unidirectional 
microscopic transitions for t — > 00. Note that if 7, = 
then H is indeed minimized by solving dH/dfa = —2ji = 
0, whereas for 7,; 7^ the minimum of H is achieved by 
taking fa as large (if ji > 0) or as small (if 7, < 0) as 
possible, i.e. fa = fa, m ax (resp. fa = fa, m in) if It > 
(resp. 7; < 0). 



Physical meaning of H 

For a start, notice that the change in yi is expressed 
in equation ^ as the sum of two terms. The first ac- 
counts for the coupling of i to other reactions in the net- 
work via shared compounds (the coupling coefficient 
is non-zero only if i and j have a metabolite with finite 
concentration in common) . Consider a species \x and two 
reactions i and j such that & > (j produces /x) and 
£f < (yu is consumed by i). Then > and a posi- 
tive net advancement of reaction j will contribute to the 
increase of yi, see Now looking at Q it is reason- 
able to expect that the probability of observing a forward 
transition for reaction i between time t and time t + 5t 
will be larger the larger (and positive) is yi{t) (in agree- 
ment with the fluctuation theorem [S]; we shall see an 
explicit example of this later on). Therefore in this case 
a positive net advancement of reaction j will ultimately 
increase the probability of a concomitant advancement 
of reaction i. In other words, this situation favors the 
emergence of a positive correlation between reactions i 
and j. Likewise, if both i and j are either producing 
(£f > 0, > 0) or consuming (£f < 0, < 0) species 
fx, their coupling will tend to anti-correlate Vi and Vj. 
In this case the dynamics discourages over-production or 
over-consumption of a chemical species. A similar pic- 
ture holds for the coefficients hi, which are related to the 
presence of sources (like nutrients) and sinks (e.g. out- 
takes) in the network. A non-zero hi acts as a force of 
magnitude proportional to that tends to polarize in 
a particular direction a reaction i that is stoichiomet- 
rically connected to a compound fi with 7^ 0. The 
favoured direction depends on whether fi is a source or 
a sink and on the sign of This effect can propagate 
to other nodes connected to i if lu^l is sufficiently large. 
The role of the concentrations appearing in and 
hi is in essence that of modulating the strength of the 
couplings and of the forcing fields with the availability of 
the intermediate metabolites. Indeed, J^ 's get stronger 
when the intermediate compounds are present in smaller 
amounts, stressing the emergence of positive correlations 
between processes (if Jij > 0) or the limits imposed by 
competition for a limited resource (if Jij < 0). When the 
concentration of the intermediate is large, instead, the 
coupling gets weaker and i and j may become effectively 
independent. The impact of hi is understood along sim- 
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thought to have 




I I >-Q O consumed by Q 

O »□ □ produced by O 

intake of □ 
Q")> outtake of Q 



iO — Oi J y >0 
iO-Oi 



Ji|<0 
h i >0 
hi <0 



FIG. 1. (a) Toy reaction network with reactions repre- 
sented as circles and chemical species as squares. Continuous, 
dashed, incoming drawn and outgoing drawn arrows denote 
stoichiometric coefficients and uptakes, respectively > 0, 



< and u M > 0. (b) Reduced reaction net- 



work with couplings and "fields" given by Q. Continuous, 
dotted, incoming grey and outgoing grey arrows denote re- 
spectively Jij > 0, Jij < 0, hi > and hi < 0. For instance, 



"£i £3 /x < 0, hi = ffu /x A > 0. Grey arrows 



are double-headed when the sign of h depends on the precise 
values of stoichiometric coefficients and uptake fluxes. For 
instance, the value of /12 = £,2 u D /x D +^2 uF / xF depends on 
the choice of the £'s and it's, since the first term in the sum 
is negative while the second is positive. 



ilar lines. In this way, the original bipartite network of 
reactions and metabolites can be re-cast as a system of 
interacting reactions with Hebbian couplings, as shown 
in Fig. [T] This is strongly reminiscent of Hopfield models 
of neural networks. In such a scenario, finding the steady 
state(s) of ([7]) is equivalent to finding the ground states 
of a system of reactions interacting with "energy" H . 

From a physical standpoint, H quantifies the resource 
mis-usage by the network so that, by minimizing H, the 
system strives to reach states in which compounds are 
used as optimally as possible, given the initial conditions 
Ui(0) (that also account for the standard GEs and hence, 
to some degree, for the a priori reversibility), the stoi- 
chiometry, the available nutrients, the production goals, 
etc. Whether for a given network the minimum of H 
is zero or not then depends on several factors, including 
the bounds on fluxes and the specific form of the uptakes. 
Note that min H = would imply that at stationarity 



(15) 



i.e. that in the steady state a Kirchhoff-law type of 
scenario holds in which strict mass-balance conditions 
are satisfied for each chemical species. The network in 
this case organizes the fluxes so that consumption and 
production exactly match for each species and meet the 
nutrient availability and outtake requirements described 
by the vector On the other hand, the physically 

relevant states with mini? > can more generally be 



V// 



(16) 



otherwise the system could e.g. consume a nutrient in 
excess of its availability. Both sets of conditions have 
been employed for the modeling of cellular metabolic net- 
works. In particular equations of the type of (15) are 



the basis of the highly successful flux-balance-analysis 
(FB A) pJH [20] , where biological functionality is included 
as an additional ad hoc constraint usually represented 
by the maximization of a specific score function (e.g. 
biomass production for E. coli in an optimal environ- 
ment). The conditions (16) are instead reminiscent of 
Von Neumann's model of reaction networks [25], where 
self-consistent flux states with a net positive production 
of intracellular metabolites can be allowed. This type of 
approach can be helpful in analyzing the metabolic capa- 
bilities of an organism and, if statistically robust produc- 
tion profiles emerge, in inferring (rather than postulat- 
ing) cellular objective functions, with the idea that chem- 
ical species that are globally produced (e.g. amino acids) 
are employed in macromolecular processes (like protein 
formation) that are not encoded by the reaction stoi- 
chiometry. We remark however that within the above 
setting finding the NESS of the system means minimiz- 
ing H , which obviously is a priori different from solving 
(fl5l) or (fl6l). 



A broader physical insight can be obtained by studying 
how H relates to standard quantities used to characterize 
non-equilibrium behavior in reaction networks, such as 
the entropy production. Indeed the entropy production 
per volume Smt of the system enclosed in the "cell" of 
volume V is given by 



TS int = -^zV 



(17) 



where x^ is the concentration of species /1 and g M = g$ + 
RT log x^ is its chemical potential. Taking time deriva- 
tives keeping in mind that at stationarity x^ is constant, 
one easily finds that H ~ ^] (i (i'') 2 /x^ — —Si nt /R. Now 

since H > 0, S,- n t < 0. In addition, H is constant in 
NESS. It follows that for time scales shorter than x^/x^ 
(i.e. for time scales over which chemical potentials are 
roughly constant) one has 



Sint — Si n t(0) — RHt 



(18) 



H is thus seen to play the role of the rate at which the 
entropy production changes in a NESS. If mini/ = 
(i.e. if the NESS is flux-balanced), then Smt must vanish 
as well (since TS mt = -£„ - «")), lead- 
ing to a state with constant entropy. If mini? > 0, in- 
stead, the entropy production is non-zero and decreases 
at the smallest allowed rate. Correspondingly, the en- 
tropy "slows down" quadratically over the time scales 
for which the theory holds. (Note that for sufficiently 
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long times the entropy production can become negative 
if H > 0, implying that this limit may be unphysical.) 
In summary, the NESS that can be obtained are cither 
characterized by H = 0, zero entropy production and 
hence constant entropy, or by H > 0, positive entropy 
production decreasing in time as slowly as possible and 
entropy increasing in time accordingly. 

It is noteworthy that this scenario essentially charac- 
terizes the Lyapunov condition described in [3D] for the 
stability of stationary states (Chapter 18). The Lya- 
punov function in our case can be computed explicitly, 
takes the form (13), its minimization bears a further 



physical meaning in terms of optimal resource allocation, 
and provides an equivalent description of the reactor as 
a system of processes interacting via Hebbian couplings. 
Indeed the quantity H appears in the Gibbs theory of 
thermodynamic stability for ideal systems. Let us con- 
sider a system at equilibrium with vanishing net fluxes 
and evaluate the effect of a perturbation that drives the 
system away from equilibrium. Such a perturbation cor- 
responds to a (small) change in the turnover Szi of each 
reaction i, which can be achieved e.g. by forcing a (small) 
non-zero flux fa through each reaction i for a time St, so 
that Szi oc faSt. The free energy per volume associated 
to the perturbed state is easily found to be given by 



G 



G eq + RTHSt 



> G, 



eq 



(19) 



in agreement with the second law of thermodynamics. 
After turning off the perturbation, the system will relax 
back to equilibrium by minimizing G (i.e. H). 



Such a transition can be fully investigated within the 
following, simple model: the stoichiometric coefficients 
are independent, identically distributed random 
numbers (such that each reaction uses and produces a 
finite fraction of the possible compounds), and the dy- 
namics of the network advances in discrete time steps of 
size St. In specific, at each time t the u^s take on values 
stochastically in {—1,1} with 



Prob{t/ i (t) = 1} 
Probst) = -1} 



a Vdt) 



n *"(*)■ 



(20) 



i.e. with Prob{^(t) = ±1} cx exp[±y l (t)/2}. Note that 
the above probability ratio is proportional to the ratio 
between the concentration of substrates and that of prod- 
ucts of the reaction. Clearly ( 20 ) does not allow to cap- 



ture the temporal structure induced by the Arrhcnius 
law, because by forcing each reaction to take place at 
every time step it neglects the fact that activation ener- 
gies (and hence characteristic timescales) can differ sig- 
nificantly across reactions. It is however reasonable to 
think that the steady state will be unaffected by these 
transients. On the other hand (20) makes the theory 



considerably easier from a mathematical viewpoint and 
the final result for the steady state (which is the focus of 
the present work) more transparent. Setting V = St = 1 
for simplicity, so that — 1 < fa < 1 for each i, one can 
follow [3T| to derive the continuous-time limit of ^ for 
large N and M. The result is the Langevin process 



= hi + 22 J i3 tanh(y i /2) 



(21) 



RESULTS 

A toy model: random reaction networks 

A simple algebraic argument allows to understand 
that, generically, a qualitative change in the solutions 



of ( 14 ) is expected to take place as one varies the ratio 
between the network parameters N and M. Let N ICV 
denote the number of reactions that remain asymptoti- 
cally bidirectional, i.e. such that jj = 0. The conditions 
Jii4>3 — that they must satisfy form a set of N rev 
equations, M of which at most are independent (strictly 
speaking, the rank of the matrix {J^} equals that of 
the stoichiometric matrix {£f }, see ([9j) ; it is however al- 
ways possible to eliminate dependencies in the latter so 
as to attain full rows rank) . Neglecting the fact that vari- 
ables are bounded, we can say that when the number of 
variables (iV rev ) exceeds that of equations the system is 
underconstrained and multiple solutions occur. So mul- 
tiple steady states exist when M < N xev . Note that iV rev 
is determined dynamically by |7| with initial conditions 
2/i(0). This suggests that different choices of j/i(0) will 
lead to different steady states, i.e. that when M < N Tev 
ergodicity (i.e. independence of the steady state on ini- 
tial conditions) won't hold. A transition is thus expected 
to occur when N Tev = M. 



where r\i is a Gaussian noise with zero mean. Clearly, 
time-averaging leads back to ( 12 ) with fa = (tanh(j/j /2)) . 
Recalling that — /3AGi = log(fa t +/fa t -) and not- 
ing that this implies fa = fa t + — fa.- = (fa.+ + 
^,_)tanh(-j9AGi/2) where, by '{20), fa\ + + fa_ 
this in turn suggests the relation 



1. 



AG, 



arctanh (tanh(j/j/2)) 



_1 1 + (tanhfa/2)) 
P ° g 1 - (tanh^/2)) 



(22) 



which explicitly links the thermodynamic driving force of 
a reaction to the (stochastic) dynamics of the quantity 
yi. The fact that the final result differs from the naive 
intuition — ftAGi — (y,) is a direct consequence of the 
stochastic fluctuations encoded in (20). Note that the 



dynamics (20) thus converges to NESS (minima of H) 
that are thermodynamically feasible, in agreement with 
the second law of thermodynamics. We can study the 



linear stability of (21) by setting 



y.i(t) = 2arctanh(0*) + X^t) 



(23) 



where {fa*} is a NESS and Xi(t) is a zero-average noise 
representing (small) perturbations to the trajectory. Re- 
actions with fa* — ±1, for which yi diverges, will be in- 
sensitive to Xi (t) . Hence it suffices to focus our attention 
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on the response of reactions for which (</>*) 2 < 1. To 
first order in Ai, fluctuations are easily found to obey the 
condition 



S« = -Jy[l-tanh 2 (y*/2)] , (25) 

where y* = 2 arctanh(0*). Now if all eigenvalues of the 
matrix Sij are positive the dynamical system will be lin- 
early stable as small perturbations occurred along the 
trajectories will die out in time. The term 1— tanh (y„/2) 
is clearly positive, hence it suffices to check that the 
smallest eigenvalue of the matrix — — X^CfCjV^-oo 
is positive. Assuming that £f are independently and 
identically distributed, the spectrum of can be com- 
puted, in the limit N — > oo with n — N/M finite, us- 
ing the results of [32]. For the smallest eigenvalue one 
finds X m in — a 2 [l — ^/n rcv ] , where n rev = N rcv /M and 
a 2 is a constant. Hence stability (and ergodicity) re- 
quires n rev < 1. When n rcv > 1, instead, the dynamics 
will be sensitive to small perturbations and, reasonably, 
its steady state will be selected by the initial conditions 
UiiO). The marginal stability condition n rcv = 1 coin- 
cides with the rough algebraic estimate of the transition 
point made above, i.e. iV rcv = M. 

To illustrate the scenario underpinned by the above 
theory, we have simulated the dynamics defined by 



Ui(t + 1) - yi{t) = E JijVj(t) + hi 



(26) 



for an ensemble of artificial reactors formed by M species 
and N reactions in order to analyze the dependence of its 
steady state(s) on n = N/M. Stoichiometric coefficients 
were chosen randomly, so that E* = with probability p 
and £i = ±1 with probabilities (1 — p)/2, independently 
on i and /i. Similarly, for the boundary metabolites we 
took v/ 1 = with probability q and — ±1 with prob- 
ability (1 — q)/2, independently on ji. For sakes of sim- 
plicity, we set — 1 for all (i (it is clear that the choice 
of xQ, does not affect the transition point where H van- 
ishes; it only changes the value of mini? in the ergodic 
phase). No prior assumption on reversibility was made, 
i.e. all microscopic transitions can initially occur in both 
directions for each process. The coefficients and hi 
are defined as in Efy and (111 except for the fact that 
Jij is re-scaled by VN to ensure that the system is well- 
behaved when N 3> 1, while Vi(t) evolves according to 



( 20 1 . For consistency, H for this case is defined as 



H 



-T- 



(27) 



Results forp = q = 1/2 are shown in Figure [2] One sees 
that initializing the dynamics from "equilibrium" con- 
ditions yi(0) — for each i one ends up in stationary 
states with H > for n < n c ~ 3 whereas H = for 



i H, unbiased i.e. 
I P 1 , unbiased i.c 
. n REV , unbiased 
H, biased i.c. 
P., biased i.c. 
,, biased i.c. 




FIG. 2. Figure 2. Average stationary values of H, frac- 
tion of asymptotically unidirectional reactions (Pi) and rela- 
tive number of asymptotically bidirectional reactions n rov = 
N lev /M versus n = N/M obtained from (261 (with no prior 



assumption on reaction reversibility for an ensemble of ran- 
dom reaction networks with NM = I0 4 constructed as de- 
scribed in the text. Averages are taken over 200 realizations 
for each value of n. Unbiased i.c. (initial conditions) refers 
to steady states of (261 for yi(0) = 0; biased i.c. instead cor- 
respond to yi(0) — 5N for each i. Note the two phases with 
H > (n < n c ) and H — (n > n c ) as predicted. The crit- 
ical point n c coincides within numerical error with the point 
where iV rcv = 1. Finally, the phase with H — is non-ergodic: 
different initial conditions lead to different NESS, character- 
ized by different values of Pi . 



n > n c . Similarly, n rcv < 1 for n < n c while n rcv > 1 
for n > n c , so that the critical point is indeed marked 
by the condition rt rov = 1. When the dynamics starts 
from yi(0) = 5N for each z, instead, the system ends up 
in a different steady state for n > n c , as signaled by the 
different value of the fraction of asymptotically unidirec- 
tional processes at stationarity, Pi = (N — N Tev )/N. In 
the phase characterized by H > 0, the steady state is 
unchanged. These results fully confirm the theoretical 
predictions derived above. We also note that (data not 
shown), if = for each fi, i.e. if there is no exchange 
with the surroundings, the corresponding networks con- 
verge to a steady state in which fa = for each i, H = 
and Pi — for each n, i.e. to chemical equilibrium. In 
other words, expectedly, boundary fluxes induce NESS. 



Red cell metabolism 

We now turn to a somewhat more realistic case in 
which nevertheless a full analytic study of the scenario 
underlying the minimization of H is possible, namely the 
standard reduced model of the hRBC metabolism, which 
includes the glycolytic and the pentose phosphate path- 
way only (21 reactions among 30 metabolites plus two 
ionic pumps, namely ATPase and NADPHase; see Ta- 
bles I and II for the numerical details and the network 
structure used here) . They operate by consuming glucose 
(GLC) in order to, respectively, maintain the osmotic 
balance through the sodium-potassium ionic pump (AT- 
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Pase), and reduce the amount of free radicals through 
glutathione reductase (NADPHase). Moreover, at a sim- 
plified level, one can think that the final products of 
these pathways are, respectively, lactate (LAC) plus an 
exchange of K + and Na + ions, and CO2. The partition- 
ing of GLC between the two pathways depends on the 
level of oxidative stress faced by the cell. Experimen- 
tal estimates based on enzyme activities range from 70% 
or more in favor of glycolysis in unstressed conditions to 
70% or more in favor of the pentose-phosphate pathway 
under oxidative stress [57]. We want to use the theory 
described here to estimate the range of variability of the 
fraction of glucose consumed by each pathway as a func- 
tion of the oxygen concentration in the environment. 

We start by assuming (reasonably) that the steady 
state of the hRBC metabolism is compatible with 
mini/ = (i.e. flux balance). A straightforward analysis 
of the emerging equations reveal that only three of the 23 
reactions are linearly independent: we choose the glucose 
uptake m glc , the flux through the Rapoport-Leubering 
shunt (or through the enzyme 2,3-DPG mutase, DPGM, 
a key step that regulates the haemoglobin's affinity with 
oxygen) 0dpgm = <^rls and the flux through the pen- 
tose phosphate pathway (or through the enzyme glucose- 
6-phosphate dehydrogenase, G6PDH) 0G6PDH = </>ppp- 
All fluxes can be written in terms of these. In particular, 
one finds 



A 9„,GLC fflPPP , 
OATPasc = iU </>RLS 



(28) 



Now the variation of the extracellular concentrations of 
GLC, LAC, K + , Na + and CO2 due to the operation of a 
single hRBC is easily seen to be given by 



^GLC = _ M GLC 



x LAC = 2u 



GLC fPPP 



i CO2 =0PPP , i Na = 3^ATPase (29) 
X K = — 20ATPaso ■ 

In turn, for the extracellular medium one has 

_ (i GLC ) 2 (i LAC ) 2 (x c ° 2 ) 2 (x Na ) 2 (± K ) 2 

— T GLC ' ™LAC ' T C02 ' T Na ' ~K 

(30) 

Minimizing this at fixed w G and i^rls one finds that 

</> ppp = 6(1 - a)u GLC - 3(1 - &)«fe LS (31) 

where a and b are defined respectively as 

_ (a; 002 )' 1 

° ~ (a;C02)-i + ( x N a )-i + (4/9)(aJC)-i + (1/9) (aM )" 1 

(32) 

( x cQ2)-i + (l/g)^ )- 1 

~ ( X C02)-1 + ( x Na)-l + (4/9)(a;K)-l + (l/9)(a;LAC)-l 

(33) 

One sees that if the concentration of CO2 is much larger 
than the others (implying a ~ 0) and </)rls — then 



^ppp ro 6 U GLC SQ fagfc the pentose phosphate pathway 
consumes roughly all of the glucose (the factor 6 is in 
agreement with the stoichiometry of carbon atoms in 
GLC and CO2). We can therefore define the fraction of 
GLC consumption through the PPP by re-scaling <fr ppp 
by 6u GLC : 



F = 1 - a - (1 - b) 



<PRLS 
2u GLC 



(34) 



From this it is also immediately clear that, in general 
conditions for CO2, F varies between 1 — a (correspond- 
ing to 0r.ls — 0) and b — a (corresponding to maxi- 
mal 4>kls — 2it GLC ). Using the typical concentration 
values for the above metabolites in the blood, namely 
x lac _ i -3 Mj x co2 _ 2 1CT 2 M, x Na ~ 0.13M and 

x K ~ 5 10" 3 M one has a ~ 0.194 and b ~ 0.625 so 
that 0.43 < F < 0.8. Despite the roughness of the net- 
work reconstruction we employed, the bounds just ob- 
tained are in remarkable agreement with experimental 
evidence. On the other hand, using the empirical esti- 
mates w GLC ~ 3 10~ 7 M/s and RLS ^ 1.4 10~ 7 M/s we 
find F ~ 0.71. 



DISCUSSION 

In this paper we have derived a variational principle for 
the NESS of chemical reaction networks, showing that, 
for time scales over which chemical potentials can be con- 
sidered constant, stationary non-equilibrium fluxes min- 
imize the function H, see (131, in which stoichiometry, 



intakes, outtakes and concentrations appear as param- 
eters. The cost function is closely related to Hopficld 
models of neural networks. This allows to rephrase such 
networks as systems of reactions interacting via Hebb-like 
rules and subject to an external forcing provided by the 
boundary fluxes. Furthermore H bears two simple physi- 
cal interpretations. First, minimizing it amounts to find- 
ing the flux organizations that keep the overall "waste" 
of chemical species to a minimum, given the boundary 
conditions to be satisfied. Second, it equals (modulo a 
constant) the time derivative of the entropy production 
per volume: NESS thus correspond to flux configurations 
such that the entropy production decreases in time at the 
smallest allowed rate, in line with the ideas exposed in 
[3D]. We have investigated the implications of such a 
picture in toy (random) chemical networks - where the 
dependence of the solutions on the network's structural 
parameters can be fully explored, revealing the existence 
of a phase transition between flux balanced states with 
min H = and states where unbalances emerge - and in 
a small real biochemical network, namely the metabolic 
network of the human red blood cell, a system that is 
possibly the closest to the theoretical situation described 
here. 

Making contact with stoichiometric models of 
metabolic networks is relatively straightforward as long 
as boundary fluxes are taken to be fixed and one does 



not include additional objective functions that NESS are 
required to maximize. For instance, the maximization of 
biomass flux (a frequent optimization criterion for bacte- 
rial metabolism }35| ) provides a source of entropy produc- 
tion even if one focuses on states with H = (as in Flux- 
Balance- Analysis, FBA). On the other hand the theory 
developed here suggests that the set of local constraints 
that describes NESS is provided by the minimization of 
H rather than by taking H = 0orH>0a priori. 

An important property that NESS should possess is 
thermodynamic feasibility, i.e. they should not contain 
infeasible loops [35] [37] . Inspired by the fluctuation the- 
orem [8], we have proposed here a simple dynamical rule 
(see (20)) that ensures that the NESS obtained as the 
minima of H are indeed void of cycles. In general (i.e. 
when straightforward minimization of H is carried out), 
it is possible to get rid of infeasible cycles by complement- 
ing the variational problem described here with the min- 
imization of the square norm of the flux vector. This is a 
consequence of the Gordan theorem of alternatives [55] : 
assuming that the matrix A = (A^) has full rank, only 
one of the following systems has a non-trivial solution: 
(a) E M 4V < Vz (for {g*} real); (b) A? h = V/i 
(for {ki > 0}). For a reaction network with stoichiomet- 
ric coefficients £f and fluxes 4n, defining A^ = <f>^ one 
sees that system (a) expresses the condition of thermody- 
namic stability ^AG, < Vi, whereas system (b) defines 
thermodynamically infeasible cycles, so that either a flux 
configuration {&} is thermodynamically feasible or it 
contains at least one cycle. Now let us consider a steady- 
state flux configuration satisfying YliZifa — uM an d let 
us assume it is thermodynamically infeasible, i.e. that 
system (b) has a solution {fc^}. We can then construct 
a new flux configuration {4>'i} as (j)^ = 4>i + Xki<f>i, with 
A a constant. Evidently, still X^ff^i — U>J ■ However 
defining Q{{(j)i}) = J2i 4>ii it is easily seen that choosing 
A so that ^Q(m) = one gets Q({#}) < Q({&}). 
In other terms, starting from a thermodynamically in- 
feasible steady-state flux configuration one can construct 
another steady-state flux configuration whose total flux 
is lower. In turn, Q has to be minimum when the flux 
configuration is thermodynamically feasible. It would be 
interesting to find a dynamical justification of this cri- 
terion. (Notice that flux minimization is an optimality 
principle frequently used for metabolic network model- 
ing, see e.g. [55].) 

Perhaps not too surprisingly, the formal aspects of this 
theory present many common traits with those devel- 
oped over the past decade for the analysis of large games 
with heterogeneous interacting agents, specifically with 
Minority Games [2H[3S]. F° r obvious reasons, such an 
analogy shouldn't be stretched and we do not elaborate 
in detail on it here. In short, however, reactions (or, 
more properly, enzymes) in chemical networks can be 
thought to be involved in a competition for the use of a 
set of possibly limited resources (the different substrates) . 
Whether a reaction can operate or not depends on how 
much substrate is available to it, i.e. on the overall sub- 



strate concentration and on how many other enzymes can 
bind the same substrates. Each reaction disposes of two 
'strategies' (or ways to access the set of resources), corre- 
sponding to the vectors of its input and output metabo- 
lites in the forward and reverse direction, respectively. 
Such strategies are anti-correlated: substrates in the for- 
ward direction are products in the reverse, and vice- versa. 



Rules like (26) can then be read as 'learning' processes 



through which enzymes try to anticipate at each time 
step whether the substrates needed for the forward or re- 
verse processes are most likely to be available, in order 
for it to operate. This parallel provides an elementary 
quantitative flavor to the idea that enzymes in biochem- 
ical reaction networks compete for the substrates. In 
absence of boundary fluxes {hi = for each i in ([T])), the 
situation is completely equivalent to the Minority Game 
with anti-correlated strategies studied in [40]. In this 
case, the dynamics converges to H = and (pi — for 
each i, i.e. the system asymptotically reaches chemical 
equilibrium. As it should be, NESS are induced by non- 
zero boundary fluxes, i.e. by non-zero fields hi. This 
additional term turns out to be the main difference be- 
tween standard anti-correlated Minority Games and the 
systems discussed here. 

Besides a possible theoretical interest in deepening the 
analogy just described (e.g. by extending the dynamical 
approaches employed for the analysis of multi-agent sys- 
tems [H] to models of chemical reaction networks) , it will 
be interesting to see how well the variational principle 
( 14 ) describes flux states in real biochemical networks. 
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NADP 


Nicotinamide adenine dinuclcot idc phosphate 


3 ± 0.5 - 10 -5 


NADPH 


Nicotinamide adenine dinuclcotidc phosphate (R) 


6±2- 10 -5 


H* 


Hydrogen ion 
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Inorganic phosphate 


1.0 ±0.5- 10" 3 


CO?, 


Carbon dioxide 
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H 2 Q* 
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solvent 



TABLE I. Metabolites appearing in the reduced model of hRBC metabolism. For reference, we include their estimated in- 
tracellular concentrations (when available; data were extracted from the BioNumbers database: see R. Milo et al. (2010) 
BioNumbers-the database of key numbers in molecular and cell biology. Nucl. Acids Res. 38 (suppl. 1) D750). The compounds 
marked with an asterisk can be subject to uptakes. 



Nr 


Abbr 


Enzyme 


reaction 


1 


HK 


Hcxokinasc 


GLC + ATP -» GQP + ADP + H 


2 


PGI 


Phosphoglucoisomcrasc 


G6P o F6P 


3 


PFK 


Phosphofructokinase 


F6P + ATP -> FDP + ADP + H 


4 


ALD 


Aldolase 


FDP *-> GA3P + DHAP 


5 


TPI 


Triose phosphate isomerase 


DHAP ^ GA3P 


6 


GAPDH 


GLyceraldehyde phosphate dhydrogenase 


GA3P + NAD + Pin 13DPG + NADH + H 


7 


PGK 


Phosphoglycerate kinase 


13DPG + ADP o 3PG + ATP 


8 


DPGM 


Diphosphoglyccromutase 


13DPG 23DPG + H 


9 


DPGasc 


Diphosphoglycerate phosphatase 


23DPG + H 2 0^ 3PG + Pi 


10 


PGM 


Phosphoglyceromutasc 


3PG <-> 2PG 


11 


EN 


Enolase 


2PG o PEP + H 2 


12 


PK 


Pyruvate kinase 


PEP + ADP + H ~¥ PYR + ATP 


13 


LDH 


Lactate dehydrogenase 


PYR + NADH + H LAC + NAD 


14 


G6PDH 


Glucosc-6-phosphate dehydrogenase 


G6P + NADP -> QPGL + NADPH + H 


15 


PGL 


6-phosphoglyconolactonasc 


QPGL + H 2 0^ 6PGC + H 


16 


PDGH 


6-phosphoglycoconatc dehydrogenase 


6PGC + NADP -> RL5P + NADPH + CO2 


17 


R5PI 


Ribose-5-phosphatc isomerase 


RL5P i-t R5P 


18 


X5P 


Xylulosc-5-phosphatc cpimerase 


RLhP o X5P 


19 


TKI 


Transketolase I 


A5P + «5P S7P + GA3P 


20 


TA 


Transaldolase 


GA3P + 57P ^ E4P + F6P 


21 


TKII 


Transketolase 


X5P + E4P o F6P + GA3P 


22 


ATPase 


Na-K pump 


ATP + H 2 -)■ ADP + Pi 


23 


NADPHasc 


Glutathione reductase 


NADPH -> NADP + H 



TABLE II. Reactions appearing in the reduced model of hRBC metabolism. Processes 1-13 belong to glycolysis, 14-21 to the 
pentose-phosphate pathway; 22 and 23 are instead the pumps. Standard reversibility assignments (based on thermodynamic 
information) are represented by the arrows. 



